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Abstract 

Synchronization induced by long-range hydrodynamic interactions is attracting attention as a 
candidate mechanism behind coordinated beating of cilia and flagella. Here we consider a minimal 
model of hydrodynamic synchronization in the low Reynolds number limit. The model consists 
of rotors, each of which assumed to be a rigid bead making a fixed trajectory under periodically 
varying driving force. By a linear analysis, we derive the necessary and sufficient conditions for 
a pair of rotors to synchronize in phase. We also derive a nondinear evolution equation for their 
phase difference, which is reduced to minimization of an effective potential. The effective potential 
is calculated for a variety of trajectory shapes and geometries (either bulk or substrated), for which 
the stable and metastable states of the system are identified. Finite size of the trajectory induces 
asymmetry of the potential, which also depends sensitively on the tilt of the trajectory. Our results 
show that flexibility of cilia or flagella is not a requisite for their synchronized motion, in contrast 
to previous expectations. We discuss the possibility to directly implement the model and verify 
our results by optically driven colloids. 
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I. INTRODUCTION 



Coordinated cyclic beating of elastic organelles such as cilia and eukaryotic flagella serve 
a multitude of functions in living organisms, ranging from motility to fluid transport and 
polarity symmetry breaking in developing embryos pj-yj. It has long been known that the 
beating cycle of cilia has a characteristic asymmetry, with two distinct parts described as 
power stroke and recovery stroke Q] , and that the cyclic pattern could lead to metachronal 
waves (of varying kinds) [5] in dense arrays of cilia [6-8]. While the necessity of this asym- 
metry for generating symmetry breaking fluid flow or propulsion could be easily understood 
from the time-reversal symmetry properties of the Stokes equation for viscous fluid flow, 
what exactly constitute the minimal conditions for synchronization and coordination be- 
ween two or more of such cyclically beating organelles is a subject of current investigation 

fll- 

There have been a number of systematic experimental studies in a variety of systems to 
)robe whether viscous hydrodynamic interaction alone can lead to synchronization as Taylor 
10| originally proposed. The experiments, which all verify the existence of hydrodynamic 
synchronization, range from studying macroscopic model flagella in highly viscous silicone 
oil (such that the low Reynolds number condition was maintained) 111 . |12J to probing the 
relative phase dynamics in pairs of beating eukaryotic flagella 13|, ll4|, to tracking colloidal 
linear oscillators using optical tweezers equipped with feedback control 15| and light driven 



asymmetrically micro-fabricated rotors 16] . Experiments on car pet s of bacteria with active 
flagella [n| and arrays of artificial magnetically actuated cilia 18-20] have revealed collective 



effects mediated by hydrodynamic interactions, such as complex flow patterns and collective 
phase shifts. 

Theoretical studies of metachronal coordination and synchronization of cyclically beating 
organelles have been performed using models and descriptions of varying levels of complex- 
ity, ranging from simple models of coupled oscillators to actuated beads and more elaborate 



elastic filament models 
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34|. While the more realistic beating elastic filament de- 



scriptions are crucial for understanding detailed properties of metachronal waves, the simpler 
actuated bead models that typically have a minimal number of degrees of freedom could be 
useful in understanding what key ingredients are needed for hydrodynamic synchronization 
to occur, and under what conditions such dynamical states could be stable. In the course the 



studies of actuated bead models, one of the questions that have been discussed in whether or 
not beads following rigid trajectories could lead to synchronization. The discussion started 
when it was shown that two rigid helices that are rotating under constant torque cannot syn- 
chronize 1 231 whi. 



e with an added a small flexibility, say to the axis of rotation, the system 
can synchronize 24(. Further studies of the actuated bead model followed the prescription 
of always having a flexible element, and somehow this was later on erroneously interpreted 



by many authors as a necessary condition for synchronization. In an earlier publication (33], 
we showed that flexibility is not a necessity, and that beads following rigid trajectories could 
lead to synchronization provided the shape of the trajectories and the beating force profile 
satisfy certain conditions. The aim of this paper is to present a thorough discussion of how 
synchronization could be achieved for rigid trajectories in a variety of cases. 

The rest of the paper is organized as follows. In Section II, we introduce the model and 
derive the coupled-oscillator equation. In Section III, generic conditions for synchronization 
are derived by linear stability analysis, and then applied to some specific trajectories and 
force profiles. In Section IV, we discuss flow properties, especially the net flow and energy 
dissipation rate in the synchronized state. In Section V, we describe the nonlinear time- 
evolution equation for the phase difference by an effective potential, which is then used 
to determine the stable and metastable stationary states for various trajectories and force 
profiles. In Section VI, the effect of flexibility is taken into account in a model of beads 
driven by moving harmonic traps. Finally in Section VII, before conclusion, we discuss the 
implications of our results to biological systems and their direct verification by optically 
driven colloids. 



II. MODEL 

A. Dynamical Equations 

We consider a pair of rotors (indexed by i = 1,2) and assume that each is a spherical bead 
of radius a that follows a fixed periodic trajectory r-j = where <pi = <fri(t) is the phase 

variable with the period 2n [see Fig. [T]. The bead is driven by an active force Fj = -Fj 
that is tangential to the orbit and is an arbitrary function of the phase. We assume that the 
two rotors are situated in parallel to each other, and that the center points of the trajectories 
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FIG. 1. A pair of rotors with the trajectory shape specified by R(0j) (i = 1,2). Each bead is 
driven by the tangential force F(<fti). The centers of the trajectories are both on the x-axis. 

are at height h from a flat substrate. The xy-plane is taken along the substrate, with the x- 
axis parallel to the line connecting the center points, and the ^-coordinate is taken vertically 
to the substrate. 

The hydrodynamic drag force acting on the i-th bead is written in the form gi = C ' 
[v(rj) — r-j]. The friction coefficient tensor £ depends on the height z of the bead from the 
substrate. However, the dependence is 0(a/z) and we neglect it by assuming a<z. Then 
the friction coefficient tensor is expressed by the friction coefficient Co = 67177a as £ = £oL 
The tangential component of the drag force is balanced by the driving force acting on each 
rotor, namely, F,i + tj • gj = 0, where tj is the tangential unit vector of the orbit given by 
tj = r£/|rj| with = drj/d0j. Substituting the expression for the drag force with r-j = r^0j 
into the force balance equation, we obtain the phase velocity as <pi — Ui + tj • v(rj)/|r^|, 
where u)i{4>i) = Fi{4>i) / Co\ r 'i\ is the intrinsic phase velocity. The reaction force — gi exerted 
by the bead on the fluid generates the flow field 



Here, G(r, is the Green function of the Stokes equation with the no-slip boundary con- 
dition at the substrate (Blake tensor). We will give its expression in the next subsection. 
On the RHS of Eq. ([1]), we assumed |r — Tj\ ^> a and retained the leading order term with 
respect to CoG(r, r^) = 0(a/\r — | ) [35|. Substituting this into the above expression for 




(1) 
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the phase velocity, we arrive at the coupled phase oscillator equation 

4>i = U)i + 



where Gy = G(rj,i\,). 



(2) 



B. Blake Tensor 



The Blake tensor G(ri,r 2 ) is given by 36] 



G^n, r 2 ) = GlM - r 2 ) - G%{v x - ri) + 2z 2 G^( ri - ri) - 2^ 2 G^( ri - ri), (3) 

where = (x,, Zj) (i = 1, 2), and fi = (x 2 , y 2 , — £ 2 ) is the point of reflection with respect 
to the substrate, and 
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where /i, v = x,y,z with summation over repeated indices assumed, are the fields of a 
Stokeslet, source doublet and a Stokeslet doublet, respectively. To 0(zf z 2 , ZxZ^), we have 



G(ri,r 2 ) = G(xi,yi,z l ,x 2 , y 2 ,z 2 

I 

3z ± z 2 



2irr] — r% 



(xi-x 2 ) 2 (xi -x 2 )(yi -1/2) -{x\-x 2 )z 2 



\ 



(xi - x 2 )(yi - y 2 ) {yi-y2? 
(Xi - X 2 )zi (2/1 - 2/2)^1 



(2/1 -2/2)^2 




\ 



/ 



(7) 



where r^- = (xj,?/j, 0) is the horizontal component of the position vector. Note that 
G(ri,r 2 ) 7^ G(r 2 ,ri) because of the cubic terms {G^ V) G^). It is convenient to decompose 
the Blake tensor into the symmetric and asymmetric part as 

/ 2x\ 2 2x122/12 X12Z12 ^ 

2xi22/i 2 22/12 2/12212 
V a;i 2 2i2 2/12^12 / 



G^=i[G(r 1 ,r i) + G(r a ,r 1 )] = i ^|- 



(8) 



Gi 2 , a = ^[G(n,r 2 )-G(r2,ri)] 3 ^ 2 
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where r i2 = (xi 2 , 2/12, z 12 ) = r 1 - r 2 and iwi 2 = 21 + z 2 . 

The friction coefficient tensor £ of the bead at height z from the substrate is given by [3 



Lu(z) = Co 
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(10) 



lbz J lbz 

up to 0(a/z). As we stated before, we assume that the relation a <C z always holds and 
neglect the correction terms. 



C. Geometric Factor 



In the following, we will assume that two rotors have the same trajectories shape and the 
force profiles. We can write each trajectory as r 4 (0) = r iO + R(0), where r i0 is the position of 
the center, and R(0) describes the shape of the trajectory. We also assume that the center 
positions are lying along the x-axis at height h from the substrate, and are separated by 
distance d a) from each other, and that their coordinates are given by 

r lo = (O,O,/0, r 20 = (d,0,h). (11) 
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We will also denote the typical size of the trajectory by b and the typical magnitude of the 
driving force by Fq: 

|R(0)|~6, F(<f))~F (12) 

Note that r£(0) = |R'(0)|t(0) where t(0) = R(0)/|R(0)| unit tangential vector of the 
trajectory. The intrinsic frequency w(0) is given by the force profile F{4>) as 

«*-A (13) 

It is useful to rewrite Eq.© 
where 

0,-) = t(&) ■ CoG(r,(0,), r,(^)) ■ t(0,) (15) 

is a dimensionless quantity of 0(ah 2 /d 3 ), and is determined solely by the geometric config- 
uration of the trajectories (i.e., shape, orientation, distance between each other, and height 
from the substrate). Hereafter we will call the geometric factor. Note that the symmetry 
relation 

#12(01,02) = #21 (02, 0l). (16) 

holds because G(ri,r 2 ) is identical to the transposed matrix G'(r 2 ,ri). 



III. LINEAR STABILITY ANALYSIS 



A. Generic Conditions for Synchronization 



Let us now examine the stability of the synchronized state by linearizing the evolution 
equation of the phase difference 5 = 0i — 2 , which reads 



01 - 02 = W(0l) - W(02) + 



W(0l)— — — - W(0 2 )- 



H 



myi, V2j 



F(0O V ^F(0 2 

Here, we used the relation ffl6|) . Setting 0i = 0(t) + 5{t), 2 = 0(t) and linearizing Eq. (TT 
with respect to 5, we obtain the linear growth rate 



(17) 



~ = «/(0) + 



H 



121 



(18) 



Integrating (}T8"]) over the period T = d(j)j(j) in the limit 5 — > 0, we obtain the cycle- 
averaged growth rate as 



[lnF(0)]'tf 12 (0,0), (19) 



~TJ 9 u;(0)[l + # 12 (0,0)] 

-~r / 1 

where the approximation is taken to the lowest order in the coupling H 12 , and T is the 
intrinsic period defined by 

20 AA 



T °=/ (20) 







The synchronized state is stable when T < 0. Equation ( TT9|) shows that a necessary 
condition for synchronization is that both the force profile F(<f>) and the geometric factor 
H 12 (<b,(j)) are not constant. However, the latter is constant only for linear and parallel 
trajectories, as we shall see below. For other trajectory shapes, the necessary condition 
for synchronization is the non-constantness of the driving force. Equation (TL9l guarantees 
that, if a force profile F(cb) makes V positive for a specific trajectory, then the force profile 
proportional to 1/F(cb) makes T negative for the same trajectory. In this sense, we can say 
that roughly half of the possible force profiles in the functional space are capable of inducing 
in-phase synchronization. We can also state that, for any given trajectory shape R(0) except 



for the linear one 



there exists a force profile F(cb) that leads to synchronization. For 



example, the force profile 



FU) = F 



i+ r diJ(H 12 (iJ,iJ)-H 12 ) 

Jo 



(21) 



where Hi 2 is the period-average of Hi 2 {<f>, <fi), makes T negative-definite and hence stabilizes 
the synchronized state. 

B. Far-Field Limit 

Let us now consider the far-field limit in which the distance between the rotors is much 
larger than the typical size of the trajectory. Also, the height from the substrate is assumed 
to be much larger or smaller than the distance: 






FIG. 2. Examples of the force profiles that act to synchronize two beads on circular trajectories 
aligned on the x-axis. (a) F{$) = F Q [l - \ sin(20)]. (b) F{4>) = F [l + \ sin (0 + f )] . 

In these limits, the Blake tensor can be approximated by the sum of isotropic (I) and dyadic 
(D) parts as 



C0G12 ^ G I (d)l + G D (d)e x e a 



(23) 



where we have used — r 2 = —de x , and the dimensionless factors Gi(d) and Gn(d) are 
given by 



Gj(d) = G D (d) = ^ 



(24) 



in the bulk geometry {h/d^> 1) and 



Gj{d) = 0, G D {d) 



9ah 2 



(25) 



in the near-substrate geometry (h/d 1). With this approximation, the geometric factor 
ffT5l) reads 



Hnfa, h) = G r {d) + G D {d) tMuifc). 



(26) 



Note that the first term is a constant and drops off from the integral ( |T9l) . Therefore, only the 
non-diagonal part of the hydrodynamic interaction controls the stability of the synchronized 
state in the far-field limit. Let us examine some specific trajectory /force profiles in this 
limit. 
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1. Circular Trajectories 



As the first example, let us consider the circular trajectory [see Fig. QJb)] 

R(0) = 6(cos 0, sin 0, 0). 

For this trajectory, we have |R'(0)| = b and t(0) = (— sin0, cos0, 0), which gives 

1 



H 



12 1 



Gd sin^ 



cos(20) + const. 



and 



T = 9r / 27r rf0[lnF(0)]'cos(20). 
J o Jo 



(27) 



(28) 



(29) 



Note that the factor cos 20 represents the second-rank tensorial nature of the hydrodynamic 
kernel. The synchronized state is linearly stable if and only if the Fourier representation of 
lnF(0) contains a negative coefficient for sin 20. Let the Fourier representation of the force 
profile be 



F(0) = F 



1 + 22 ^ n Sm ( n + $n 



n=l 



(30) 



where we assume < A n < 1 to avoid singularity of InF. Up to O(A^), we obtain the 
growth rate as 

io L n=l 

The only harmonic mode that contributes to T at 0(A n ) is n = 2, for which T is most 
negative at <5 2 = 7T. In this sense, the force profile that is most efficient in inducing synchro- 



2A 2 cos 6 2 + A\ sin(2«Ji) - ^ 2A n+2 A n sin (5 n+2 - 5 n ) 



(31) 



nization is given by 



F(0) = F [l - A 2 sin(20)], < A 2 < 1. 



(32) 



Another harmonic mode that contributes to T by itself is n = 1, for which T is most negative 
at S\ = 7r/4, 57r/4. Thus we obtain the candidate force profile 



F(0) 



l + Ai sin 



-1 < A < 1. 



(33) 



These two force profiles are illustrated in Fig. [2j Note that higher harmonic modes (n > 3) 
can stabilize synchronization only when they are mixed with the other modes. 
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Next, we consider rotated circular trajectories. The trajectory (1271) rotated around the 
y- axis by angle a, R(0) = 6(cos a cos 4>, sin</>, sin a cos0), gives the additional factor cos 2 a 
to the linear growth rate via Eq. f|26|) . Note that if the trajectory planes are perpendicular 
to the x-axis (a = vr/2), we have T = and the synchronized state is only marginally 
stable. On the other hand, rotation around the x-axis does not change the linear growth 
rate, because the hydrodynamic kernel (1231) is invariant for the rotation. However, near-field 
corrections will introduce an important dependence, as we shall see in the next section. We 
do not consider rotation around the z-axis, which is equivalent to shift of the phase by a 
constant. 

2. Linear Trajectories 
The linear trajectory 



gives t(0) = sgn[R'((j))]e x , which makes the geometric factor (1261) constant. Thus, at the 
level of linear stability analysis, the synchronized state is neither stabilized nor destabilized 
for any force profile. However, nonlinear stability analysis shows that the stability is weakly 
affected by force modulation, as we shall see in the next section. 

3. Elliptic Trajectories 
For the elliptic trajectory 



the x-component of the tangential vector t x ((f>) = b x cos 0/ y b\ cos 2 <p + 6 2 sin 2 contains all 
the harmonic modes with odd n if b x ^ b y . As a result, force modulations containing any 
harmonic mode with even n can induce synchronization at 0(A n ), if the Fourier coefficients 
are suitably chosen. For example, when b x > b y , the force profile 



R(0) = R{<f>)e, 



(34) 



R(0) = (b x cos(j), b y siii(f)), b x ,b y >0, 



(35) 




F{(j>) = F [1 + A 4 sin(40)] , < A 4 < 1 



(36) 




X = 




(37) 
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Force modulations with odd harmonics can also induce synchronization, because they give 
rise to even harmonic modes in lnF(0), but only at O(A^). 



IV. FLOW RATE AND ENERGY DISSIPATION 



Now let us see how synchronization affects flow properties in the substrated geometry. 
We define the volume flow rate Q as the flux through a half-plane in the "down stream" 
(x — > — oo): 

/OO POO 
dy / dzv x (r,t) 
-oo Jo 

/OO POO 
dy / dz G xv {r,Yi) ■ g { . 
-oo Jo , 



(38) 



In the second line we used the expression for the flow field (pQ). Note that, due to volume con- 
servation, the integral in (|38|) does not depend on the x-position of the half-plane. However, 
it is easier to calculate it in the limit x — > — oo, where we can use the 0(r~ 2 ) approximation 
for the Blake tensor, 



Gfiu{ T , r i) 



3zzi 



xy y 2 yz 




(39) 



which gives 



Q(t) = —Y][h + R g {(f> i )]g ix . 



(40) 



Because we are interested in the change in the flow rate due to hydrodynamic interaction 
between the rotors, we retain the first order term with respect to G in calculating the drag 
force, which reads, 



Si = Co[v(r. ; ) - ii 
= Co 



(41) 



CoG(r,, rj ) ■ R%)4- - R'(0i)0i 

In the in-phase synchronized state fa = <pj = <p, the two rotors have the same period T, and 
the cycle-averaged flow rate is calculated using ( HOI and (!4TI) as 



Q = ^ 



T 



dtQ(t) 



Qa 



E 



2tt 



[h + RMl 



(42) 
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We shall use the far-field approximation (|23|) and fl25|) . to obtain 



12a(l-G D . 



2- 



Q = rp 1 I d^R z {(j))B! x {4>). (43) 



o 



Note that the flow rate is zero for planar geometry (R z ((f)) = 0). The hydrodynamic inter- 
action modifies the flow rate not only through the prefactor 1 — Go but also through the 
period T, which is given by 



2i i / r 2iT 



T= -f- -±si l ~ H MM- (44) 

Jo Jo ^W) 
We can also calculate the power needed to drive the beads, which is given by 

P(t) = X>.(- ft ). (45) 

i 

Its cycle-average in the synchronized state is calculated to the first order of G as 



2tt 



T 
T Jo 



o 



|R(0)| 2 - u;(0)R(0) • ^CoG(r,, ri ) • R(0) 

r dxfruiW {[1 + H l2 {<t>,4>)\ |R(0)| 2 - G D R'M 2 } ■ (46) 
Jo 

For example, let us compute the flow rate and power for the vertical circular trajectory 

RO) = 6(cos0,O,sin0). (47) 

For this trajectory, the integrals in (T4"3"j) and f !4B|) give nb 2 and 2nF b/( , respectively, where 
F Q is the cycle-average of the driving force. It yields the mean flow rate 

Z}= iwq -Go) m 

and the mean power 

P = -y 2 -- (49) 
The period T depends on the force profile, and is given by 

T= 2 -^(t -t iGd ), (50) 

-TO 

1 f 2w „ Po 



ro= ^y (51) 

ri ^y (52) 
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FIG. 3. The coefficients To and n in Eq. ([52|) for the force profiles ([32]) and ([33]) . as functions of 
A% and (resp.). For the two force profiles, the curves are identical and t\/tq is equal to 1/2. 

The dimensionless coefficients Tq and T\ are positive for any force profile (with F(<f>) > 0), 
which means that the period decreases by the hydrodynamic interaction. Furthermore, we 
have ti/tq < 1 for any force profile, which means that the mean flow rate also decreases by 
the hydrodynamic interaction. 

In Fig. El we plot To and t% for the force profiles ( 1321) and ( 1331) as functions of the 
amplitude A 2 and A 1: respectively. The dependencies on A\ and A 2 turn out to be identical 
for each of r and 7~i, with the ratio t\/tq equal to 0.5. Both of the coefficients, and hence 
the period, diverge as A\^ are taken to unity. When A\^ = 1, there are stall points (where 
F((j)) = 0) on the trajectory, and it takes infinite time for the bead to pass the points. 

V. NONLINEAR ANALYSIS 

In this section, we analyze the fully nonlinear evolution equation for the phase difference 
(JTTj), which allows us to explore various dynamical states and their stability. By using the 
full Blake tensor, we will also discuss the near-field effects due to finite size and height of 
the trajectories. 

A. Effective Potential 

The difference in the phase velocities (TT7T) consists of the intrinsic phase velocities (the 
first and second terms in the RHS) and the interaction term (the third term in the RHS). In 
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order to focus on the latter, we exploit the gauge invariance of Eq. flTT]) . i.e. the invariance 
under the transformation — > $(0) where $ is a new phase variable (or a "gauge") satisfying 
$(0 + 27r) = $(0) +27T. We choose the specific gauge $ that gives a constant intrinsic phase 
velocity, which we will call the canonical gauge. It satisfies $ = 2tt/T = Q in the absence 
of hydrodynamic interaction, and is obtained from the original gauge <fi via the relation 

.ColR'MI 



— - - - <r 



F(<j>) 



(53) 



In the canonical gauge, the intrinsic terms in Eq. (I17j) cancel out, and the phase difference 
A = $! - $ 2 obeys 



A = Q 



F($ 2 ) F($i 



#i 2 ($i,$ 2 ), 



(54) 



where the force profile F(Q) and the geometric factor if 12 ($i, $ 2 ) are related to those in the 
original gauge via F(<3>) = F{4>) and #i 2 (<I>i, $ 2 ) = H 12 ((pi } 2 ). Note also that R($) = R(0) 
and R'(<I>) = ^|R'(0) = ^-F(0)t(0). We rewrite ( 1541) in terms of A and the phase sum 
S = $ x + $ 2 , as 



n 



#12 ( 



S+A S-A 
2 ' 2 



) =W(E,A). 



(55) 



Note that A/f2 = 0(Gr>) <C 1, while E/f2 is an 0(1) quantity. Therefore, we can approxi- 
mate A to be constant over one period where £ increases by 4it. With the approximation, 
we take the average of (|55|) over one period < t < T, to obtain 

17T 



A 



1 

4?r J() 



ciS W(E, A) = W(A), 



(56) 



which defines the effective force W(A). We introduce the effective potential V(A) by 

■A 



V(A) 



dA'W(A') 



(57) 



with which the dynamics reduces to minimizing the potential: 

A - dV 



(58) 



Thus we have eliminated the fast variable £ and describe the slow dynamics only by A. 
This approximation is shown to be to correct to the lowest order in the interaction 
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Note that the factor F (^) /F p±^) -F (^)/F (^) in Eq. (ESJ is an odd function 
of A by construction. The other factor H^^^-^^— ) is an even function of A if the 
following identity holds: 

#1 2 ($1,$2) =#l 2 ($2,$l). (59) 

This is the case for the far-field limit (6/d — > 0), as we see from Eqs. (ll5l) and (1231) . In that 
case, W^(£, A) and V^(A) are odd functions of A and hence the effective potential is an even 
function: V(A) = V(-A). 



B. Far-Field Limit 



First, let us derive the effective potential in the far-field limit (b/d — > 0), for trajectories 
in the bulk (h/d — > oo) and/or near substrate (h/d <C 1). 



1. Circular Trajectories 

First let us consider circular trajectory (J27J) with the force profile (1321 . The phase in the 
canonical gauge is obtained via fl53|) . as 

m *W*r^w < 60 > 

K{27r) J 1 — sin 20' 

Accordingly, the intrinsic phase velocity in the canonical gauge is given by 

fi= K(27)' fi ° = ^ (61) 

For A 2 < 1, we can approximate 0($) and -F($) as = $ + 4f cos 2$ and -F($) = F($) to 
0(yl 2 ), which gives 

V(A) ~ V (A) = n G D A 2 (l - cos A). (62) 

For not small values of A 2 , we compute the integrals in (1561) and (I6"0"j) numerically. We plot 
V(A) in FigJH which shows that the approximation ([6"2l is very good. Even for A 2 = 0.99, 
the deviation |V(A) — Vo(A)]/Vo(A) falls within 11 % for any value of A. Also note that 
V(A) is an odd function of A 2 . For A 2 > 0, it is minimized at the in-phase synchronized state 
(A = 0), while for A 2 < 0, it is minimized at the anti-phase synchronized state (A = ir). 
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FIG. 4. The effective potential V(A) in the far-field limit b/d — > for the circular trajectory (|27p . 
either in the bulk (h/d ^> 1) or near the substrate (h/d <C 1), for (a) the force profile (|32p . and 
(b) the force profile 

For the force profile (|33|) . the effective potential can be calculated in a similar way as 
above, and is plotted in Fig. HI A perturbative calculation to 0(A\) gives 

V(A) ~ ^o(A) = ^° G 4 ° A ' (1 - cos A). (63) 

Again, the approximation ( 1631) is good for moderate values of Ai. The deviation [V(A) — 
Vo(A)]/Vo(A) falls within 0.02 for A\ = 0.5. Although the deviation increases to 0.57 at 
A\ = 0.99, the shapes of V(A) and Vq(A) are quite similar. 



2. Linear Trajectories 

Next let us consider the linear trajectory fl34|) with R(<f>) = b cos <fi using the near-substrate 
approximation. In Fig. [51 we plot the effective potential for (a) the force profile ( 132]) and 
(b) the force profile fl33|) . In both cases the potential is minimized at A = 0. This result is 
not expected from the linear stability analysis, which showed that the in-phase synchronized 
state is only marginally stable for any force profile. It suggests that the potential scales with 
V(A) oc |A| K with the exponent k > 2 near A = 0. In FigJBJc), we plot V'(A) = -W(A) 
for the case (a). It indicates the non- analytical behavior k = 5/2, which we will prove in 
the following paragraph. 
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FIG. 5. The effective potential V(A) in the far-field limit b/d — > for the linear trajectory ([2 
near the substrate (h/d -C 1), for (a) the force profile ([32]) and (b) the force profile ([33]) . (c) 
Non-analytic behavior V'(A) = — W(A) oc A 3 / 2 for the case (a). 

The gauge condition fl53|) gives 
We see that $(77 + (f>) = n + $(0), and especially $(-7r) = 7r. For |0| <C 1, we have 
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K{4>) ps sgn(0) • 2 /2 and hence $(0) ps [7r/iT(27r)] sgn(0) • 2 . 
We also expand the factor in fl55|) in powers of A as 



while the other factor behaves like a step function 
'E + A S-A 



In F I - 



+ 0(A 3 



(65) 



12 



G D sgn 



sin 



E + A 



sin 



A 



(66) 



For small and positive value of A, the latter equals — Gd when Inn — A < E < 2n7r + A (n: 
integer) and equals Gd otherwise. These give the effective force to 0(A 2 ) as 



2tt 



f(-f) 



2^ 



IAI • A, 



7T 



(67) 



where we used lnF(<3>) = lnF(0) ps lnF - 2A0 p^ lnF - sgn($) • Ay/[K(2iv)/ir]\®\, which 
is an approximation for |<E>| 1. Thus we obtained the non-analytic behavior W(A) oc 
-sgn(A) • |A| 3 / 2 , or V(A) oc sgn(A) • |A| 5 / 2 . 



3. Elliptic Trajectories 

Next we consider the elliptic trajectory ( I35jl in the near-substrate approximation. For 
elliptic trajectories, the tangential vector t(0) = R'(0)/|R'(0)| and hence the geometric 
factor fT26|) contain various harmonic modes, which produce richer behaviors than the circular 
trajectories. In Fig. [61 we show the effective potential for (a,b) the force profile ( |32|) with 
A 2 = 0.5, (c,d) the force profile ( l33l) with A\ = 0.5, and (e,f) the force profile ( 1361) with 
Ai = 0.5. In (a,c,e), we show the potential curves for b x > b y , while in (b,d,f), the potential 
curves for b x < b y are scaled by (b x /b y ) 2 . (Note that the potential converges to zero in the 
limit b x /b y -)■ 0.) 

In (a), the potential has a single minimum at A = 0. As b x /b y — > 0, the scaled potential 
V(A)/(b x /b y ) 2 converges to a V-shape curve. 

In (b), a local minimum at A = ±7r appears for b x /b y < 1 in addition to the minimum 
at A = 0. For b x /b y ^ 0.6, the anti-phase synchronized state becomes stable, while the 
in-phase state becomes metastable. For 62 = 37r/4, the sign of the potential has an opposite 
sign, and we obtain bistable minima at A = ±A , with A — 7r/2 for b x /b y = 0.4 and 
A — )• 7r as b x /b y is increased to unity. 
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FIG. 6. The effective potential ^(A) in the far-field limit b/d — > for the elliptic trajectory (|35p 
near the substrate (hjd <C 1), for (a,b) the force profile ([32]) with A2 = 0.5, (c,d) the force profile 
(|33|) with A\ = 0.5, and (e,f) the force profile (|36|) with = 0.5. The long-axis of the ellipse is 
along the x-direction in (a,c,e), and along the y-direction in (b,d,f). 
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FIG. 7. (a) The phase function $(<^>) for the elliptic trajectory with b y /b x = 0.5 and the force 
profile (|36p with = 0.5. (b) Temporal oscillation of 5 in the phase-locked state A = Ao = 0.992. 

In (c), we have bistable minima at A ~ ±A with 7r/2 < A < it for b x /b y < 1, and 
< A < 7r/2 for b x /b y > 1. A metastable minimum is located at A = and A = ±7r, for 
b x /by < 1 and b x /b y > 1, respectively. For b x /b y = 1, no phase locking occurs because the 
potential is constant. In (b) and (c) also, the scaled potential V(A)/(b x /b y ) 2 converges to a 
master curve with sharp peaks and valleys in the limit b x /b y — > (not shown). 

Now let us consider the meaning of the minimum at non-zero A. In the phase-locked 
state, the phase difference A in the canonical gauge is constant, but it generally means an 
oscillation of the phase difference 5 in the original gauge, because it is a function of both A 
and £ ~ 2Qt. Let us take for example, the elliptic trajectory with b y /b x = 0.5 and the force 
profile (|36D with A A = 0.5. 

In Fig. El^a) we show the phase function $(0), which has the period 7r. The figure also 
illustrates the relation between 5 = (f>i — <p2 and A = $x — $ 2 - The effective potential, shown 
in Figj6]^e), has double minima at A = ±Ao with Ao = 0.992. 

Also, it has a metastable minimum at A = n. In Fig. El^b), we show 5 in the phase-locked 
state A = A as a function of time (using the relation Qt = E/2, where the origin of time 
is chosen arbitrarily). It oscillates with the period of $(0). The amplitude of oscillation is 
as large as 1.25, and is larger for stronger modulation of the force profile (that is, for larger 
amplitude A±). On the other hand, 5 remains constant in the anti-phase synchronized state 
A = 7r, because it matches the period of the phase function. 
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C. Near-Field Corrections 



Next we consider the near-field effects arising from finite size of the trajectory, by using the 
full Blake tensor G12 given by The finite size of the trajectory introduces dependences 
of G12 = Gi 2 (r 10 + R($ 1 ), r 2 o + R($2)) on $! and $ 2 - in general, Gi 2 is not symmetric with 
respect to the exchange of $1 and $2, which leads to asymmetry of the effective potential 
V(A). To see this explicitly, we expand the Blake tensor to the first order with respect to 



R(gj 
d 



1,2, 



(68) 



which is assumed to be small. We also assume that the height of the trajectory is of the 
same order as its size, and introduce the dimensionless height, 

// = '-. (69) 



h 

7r 



Then we can use the 0(h 3 ) approximation (J7J) as a starting point for the expansion. Sub- 
stituting X12 = d(-l + p lx - p 2x ), 2/12 = d(pi y - piy), Zi = d(h + p iz ), Z12 = d(p lz - p 2z ), and 
W12 = d(2h + pi z + P22) into (!H|9|) and retaining 0(p, h) terms, we have 
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When we exchange 0i and 02 or p l and p 2 , the sign of Gi2, s is reversed, while Gi2, a remains 
unchanged to this order. 

In the bulk geometry h/d — > 00, the effects of finite trajectory size can be examined more 
simply, by expanding the Oseen tensor to the first order in p. It gives 

' 2(1 + p lx - p 2x ) -(ply ~ p2y) ~{plz ~ P2z) ^ 

\Ply-P2y) l+Plx~P2x (72) 

\Plz-P2z) l+Plx~P2z 



CoG(n,r 2 )~C, 



hulk 



\ 



J 
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with Cbuik = 3a/ Ad. Therefore all the terms change their signs upon exchanging <f>\ and 02- 
In the following, the effective potential is calculated for the circular trajectory f[2"Tl) rotated 
around the x-axis by the angle 0, or 

R(0) = 6(cos 0, cos (3 sin 0, — sin (3 cos 0) . (73) 

The aspect ratio b/d will be fixed to 0.05 unless otherwise stated. We use both the full 
Blake tensor fl3]) and its first-order approximation [f lTUj) . flTTj) . and ([72]) ]. and compare with 
the zeroth-order (far-field) results. 

1. Force profile 132^ with A 2 = 0.5 

The effective potential is plotted in Fig. [8]for trajectories lying (a) in the bulk (h/d — > oo), 
(b) in a plane horizontal to the substrate (h/d = 0.1, /3 = 0), and (c) vertical to the substrate 
(h/d = 0.1, /3 = 7r/2). In (a) and (b), the potential is an even function of A and is well 
approximated by the zeroth order result (163]) . while in (c), V(A) has a negative average 
gradient with V(tt) < V(—n). 

2. Force profile with Ai = 0.5 

The effective potential is plotted in Fig. |9] The trajectories are lying either (a) in the 
bulk (h/d — > oo), (b) in a plane horizontal to the substrate (h/d = 0.1, /3 = 0), or (c) in a 
plane vertical to the substrate (h/d = 0.1, /3 = vr/2). In (a) and (b), the potential curves 
have negative average gradient with V(n) < V(—ir), and have local minima at A = 0. Note 
that the asymmetry of the potential curve is larger in the bulk case. In (c), the potential 
has a positive average gradient with V(tt) > V(—tt), and has two metastable minima. There 
is a saddle point at A = 0. Note that all these features are already seen in the first-order 
approximation. 

In Fig. [TUT a). we plot the potential for different trajectory size in the bulk geometry. 
The average gradient of the potential is enhanced with h/d. In Fig. [TWb). we plot the 
potential for the vertical trajectory ((3 = it/ 2) and with different height. Note that the 
average gradient of the potential changes its sign from positive to negative at intermediate 
height. In Fig. [TOT c). we show the dependence on the tilt angle (3. The average gradient 
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FIG. 8. The effective potential V(A) for the circular trajectory (I27D with b/d = 0.05 and the force 
profile (|32p with Ai = 0.5. The trajectories are either (a) in the bulk (h/d — > oo), (b) horizontal 
to the substrate (h/d = 0.1, j3 = 0), or (c) vertical to the substrate (h/d = 0.1, j3 = vr/2). Shown 
are results from the full Blake tensor as well as its zeroth-order and first-order approximations in 
terms of b/d. 
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FIG. 9. The effective potential V(A) for the circular trajectory (I27D with b/d = 0.05 and the force 
profile (|33p with A\ = 0.5. The trajectories are either (a) in the bulk (h/d — > oo), (b) horizontal 
to the substrate (h/d = 0.1, j3 = 0), or (c) vertical to the substrate (h/d = 0.1, j3 = vr/2). Shown 
are results from the full Blake tensor as well as its zeroth-order and first-order approximations in 
terms of b/d. 
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FIG. 10. The effective potential V(A) for the circular trajectory (|27p with b/d = 0.05 and the force 
profile (j33|) with A\ = 0.5. Shown are dependencies on (a) b/d in the bulk geometry (h/d — > oo), 
(b) h/d, and (c) the tilt angle /3 of the trajectories. 



of the potential changes its sign around /3 = 7r/6. We thus find that the asymmetry of the 
potential sensitively depends on the size, the height, and the tilt of the trajectories. 
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VI. EFFECT OF FLEXIBILITY 



Hitherto we have only considered rotors with rigid trajectories, but in a real system the 
trajectory could be affected by hydrodynamic flow due to flexibility or compliance of the 
rotor. As an example, let us consider a bead driven by optical tweezers, whose focus moves 
along a prescribed trajectory By controlling the distance between the focal point and the 
bead, one can tune the tangential driving force 40|. We approximate the potential created 
by the laser beam by the harmonic potential U(S) = |>S 2 , where S is the displacement of 
the bead from the focal point. The bead of the z-th rotor is thus positioned at r$(t) = 
r-jo + R(0j(t)) + Sj(i), and its velocity is r, = R'(0j)0j + Sj. The traction force due to the 
laser beam is balanced with the viscous drag force as 

feS i = Si = Co[v(r i )-r i ], (74) 

while its tangential component is prescribed as F{(f>i) = —kSi ■ t (</>;), or 

Ft = -& ■ U, (75) 

where abbreviations Fi = F(<f)i) and tj = t(0j) are used as before. In the limit k — > oo, we 
restore the model of rigid rotors developed in the previous sections. Let us derive the phase 
evolution equation by expansion in powers of k~ l . 

First, the intrinsic phase velocity cjj = u)((j>i) is determined by setting v(r) = 0. Equation 
(1741) gives 

= _ tI (r ^ + ^ ~ "tl^ - ~V F ^ w 

up to 0(A; _1 ). Using this in f J75|) . we obtain the intrinsic frequency to 0(A; _1 ) as 

(77) 

The hydrodynamic interaction is incorporated by substituting 

v(r 4 ) = -J2 • Si - - E CoG ^ ' { R '^ + ^i) ( 78 ) 

3 3 

into f!74|) . For simplicity let us assume the far- field limit, where is given by the constant 
symmetric tensor (1231) . After some calculation, we obtain the phase evolution equation in 
the form 
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where the function H 



■i j 



H, 



and Jj n — F 



Hij — tj ■ CoGjj 
ij, 4>j) is defined by 



now includes an 0(l/k) correction as 



t' 



fclR'J 



(80) 
(81) 

The stability of the synchronized state is examined by setting 0i = + <5 and 02 = and 
linearizing the evolution equation of S as before. After some straightforward calculation, we 
obtain the cycle-averaged growth rate of the phase difference 5 as 

1 r 27T 

T = — / d(P{-2[ln(\R'\cu)]'H 12 + (\ncu)'J 12 + AH[ 2 + AJ[ 2 }, (82) 
-*0 Jo 

where the functions in the integrand are to be evaluated at 0i = 02 : 



and we define 
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(83) 



T 



for any two- variable function A(<p\, 02). 

For example, let us consider the circular trajectory R(0) = 6(cos 0, sin 0, 0) = 6n(0). 
Using |R'(0)| = b, n'(0) = t(0), t'(0) = — n(0), and the Fourier representation (130]) of the 
force profile, we obtain the growth rate up to 0(A n ) as 

A 2 G D (cos S 2 - ~ sin 5 2 ■ ^ - (4G 7 + 2G D ) ^ . (84) 

This should be compared to the result (!3~T!) for the rigid rotors. We see that the flexibility 
tends to enhance synchronization due to the last term on the RHS. Note also that the 
small parameter representing the flexibility is Fo/kb. If the displacement from the focal 
point (which has the typical magnitude Sq ~ F /k) is much smaller than the size of the 
trajectory, which is the case in the optical tweezer experiment [4^], the flexibility has only a 
weak effect in inducin g sy nchronization. These results qualitatively agree with the findings 



of the previous study 29( that assumed constant driving force and radial displacement from 
a circular trajectory. In the paper, the model parameters are estimated for cilia, which give 
the dimensionless coupling (that corresponds to our Fo/kb) to be on the order of 10~ 2 — 10~ 3 . 
On the other hand, we can expect 0(1) modulation of the driving force from the effective and 
recovery strokes of cilia. Therefore, we conjecture that the force modulation plays dominant 
roles in establishing the coordinated ciliary beating. Finally, we mention that our model of 
flexibility can be also modified for rotors allowing tangential displacement, such as a bead 
attached to the tip of an elastic rod. 
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VII. CONCLUDING REMARKS 



By linear and nonlinear analysis of the coupled oscillator equation, we have fully charac- 
terized the dynamical states of a pair of rotors making rigid trajectories. In particular, we 
obtained the necessary and sufficient conditions for in-phase synchronization, which show 
that a wide variety of beating patterns induce synchronization for an arbitrary trajectory 
shape. Even for parallel linear trajectories, which predict only marginal stability in the lin- 
ear analysis, the effective potential has a global minimum at the in-phase state if we choose 
a suitable force profile. 

The results confirm and strengthen our previous finding 33| that flexibility of the rotors is 



not a requisite for synchronization, although it has been highlighted in many other studies. 
In the present paper, we incorporated flexibility into our model and explicitly compared 
its effect to the effect of force modulation. If the disturbance of the trajectory due to 
hydrodynamic interaction is small compared to the size of trajectory, the flexibility has 
only weak effect in establishing synchronization. For cilia, sizable modulation of the driving 
force is expected from their effective and recovery strokes, and it should play a dominant 
role in coordinating their beating. Recently, another mechanism for driving synchronization 



a swimming Chlamydomonas has been proposed using a simple three- 



between flagella o 
sphere model (42, 



431 ] . In these studies, the phases of the two flagella are predominantly 
coupled via translation and rotation of the cell body. This type of coupling originates 
from the condition that the net force and torque acting on the cell vanish, and is specific 
to rotors attached to a freely-suspended body. The condition for synchronization with 
this type of coupling is different from hydrodynamic synchronization (for example even 
constant forcing profile and circular trajectory could lead to synchronization under certain 
circumstances). Also, the coupling is weaker than the hydrodynamic one if the cell body 
is much larger than the distance between neighboring flagella or cilia, which is the case in 
densely flagellated/ciliated cells such as Volvox and Paramecium. 

The effective potential that governs the nonlinear dynamics of the rotors have a number 
of remarkable features. First, it allows us to locate all the stable and metastable states of the 
system at a glance. In the far-field limit, the potential is symmetric. For circular trajectories 
with simple force modulation (consisting of a single harmonic mode), the potential has only 
one minimum that describes either in-phase or anti-phase synchronization. Bistable and 
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metastable states appear for more complex trajectory shapes such as ellipses. When the 
system is trapped in an out-of-phase stable/metastable state, the phase difference (in a 
natural gauge) oscillates as a function of time. We have incorporated near-field corrections 
due to finite trajectory size, and found that the overall shape of the potential, especially 
its average gradient, sensitively depends on the size/height/tilt of the trajectory. When the 
potential has a non- vanishing average gradient, each of its local minimum corresponds to a 
metastable state. In the presence of strong noise, we may observe phase slippage in a specific 
direction. We note that experiments on the flagellar beating of a mutant Chlamydomonas 
have recently shown anti-phase synchronization [4l|. It will be interesting to probe the 
differences between this mutant strain and the wild-type Chlamydomonas in terms of the 
beating pattern of the flagella, and examine whether the phenomenon can be quantified 
within the framework of our model. 



A more direct experimental test of our findings could be pursued in a simpler system 
that does not have the complexity of the living organisms, such as optically driven colloids. 
Optical tweezers with moving focus can drive the colloidal particle on a prescribed trajectory, 
and by controlling the distance between the focal point and the bead, one can also prescribe 
the force profile. Experiments are currently underway in the group of Pietro Cicuta at the 
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have been used 



Cavendish Laboratory along these lines |40|]. Also, optical vortices 
to trap colloidal particles on a ring and drive them in one direction. The driving torque 
could be modulated by tailoring the helical structure of the laser beam to give a prescribed 
force profile. 



In forthcoming papers, we plan to discuss the collective dynamics of arrayed rotors, and in 
particular the formation of traveling waves. Such a study should become a first step towards 
understanding the relation between the beating pattern of cilia and the metachronal waves 
they form. We will also consider a pair of rotors with different intrinsic frequencies, which 
will induce phase slips similar to those observed in Chlamydomonas 13j, Il4| . The present 
paper assumes spherical beads, but the analysis could be extended to non-spherical bodies 
such as rods or helices, which are closer to the shapes of biological filaments. 
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